#include "moremath.h"
using namespace std;

namespace moremath{

vector<double> solveQuadratic(double a, double b, double c){
	vector<double> ret;
	ret.clear();

	if(a==0 && b!=0) ret.push_back(-(double)c/b);
	else if(a!=0) {
		double det=b*b - 4*a*c;
		if(det>0){
			double t1, t2;
			t1= (-b+sqrt(det))/(2*a);
			t2= (-b-sqrt(det))/(2*a);
			ret.push_back(t1);
			ret.push_back(t2);
		}
		else if(det==0){
			double t;
			t = -(double)b/(2*a);
			ret.push_back(t);
		}
	}
	return ret;
}

vector<double> solveCubic(double a, double b, double c, double d){
	vector<double> ret;
	ret.clear();
	if(a==0) return solveQuadratic(b, c, d);

	//eqn is z(3) + a2*z(2) + a1*z + a0 = 0
	double a2= (double)b/a;
	double a1= (double)c/a;
	double a0= (double)d/a;
	
	// after substituting z as x-a2/3
	// eqn becomes x(3) + px - q = 0
	double p= (3*a1- a2*a2)/3;
	double q= (9*a1*a2- 27*a0- 2*pow(a2,3))/27;

	// now substitute x as 2*sqrt(abs(p)/3)*y
	// eqn becomes 4*y(3) + 3*sgn(p)*y = C
	
	double C = (q/2)*pow(3/abs(p), (double)3/2);
	
	if(p>0) {
		double y= sinh(((double)1/3)*asinh(C));
		double x = 2*sqrt(abs(p)/3)*y;
		double z = x - a2/3;

		ret.push_back(z);
	}
	else if(p<0){
		if(C>=1){
			double y= cosh(((double)1/3)*acosh(C));
			double x = 2*sqrt(abs(p)/3)*y;
			double z = x - a2/3;

			ret.push_back(z);
		}
		else if(C<=-1){
			double y= -cosh(((double)1/3)*acosh(abs(C)));
			double x = 2*sqrt(abs(p)/3)*y;
			double z = x - a2/3;

			ret.push_back(z);
		}
		else {
			/* three solutions  ??? */
			double y1 = cos(((double)1/3)*acos(C));
			double y2 = cos(((double)1/3)*(2*M_PI + acos(C)));
			double y3 = cos(((double)1/3)*(4*M_PI + acos(C)));
			double x1 = 2*sqrt(abs(p)/3)*y1;
			double z1 = x1 - a2/3;
			double x2 = 2*sqrt(abs(p)/3)*y2;
			double z2 = x2 - a2/3;
			double x3 = 2*sqrt(abs(p)/3)*y3;
			double z3 = x3 - a2/3;

			ret.push_back(z1);
			ret.push_back(z2);
			ret.push_back(z3);
		}
	}
	else {
		if(q>0) {
			double x = pow(q, (double)1/3);
			double z = x - a2/3;
			ret.push_back(z);
		}
		else {
			double x = -pow(abs(q), (double)1/3);
			double z = x - a2/3;
			ret.push_back(z);
		}
	}
	return ret;
}

}
